{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convex approximate dynamic programming\n",
    "\n",
    "We consider a stochastic control problem of the form\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & \\underset{T \\to \\infty}\\lim {\\mathbb E} \\left[\\frac{1}{T} \\sum_{t=0}^{T-1} \\|{x_t}\\|_2^2 + \\|{\\phi(x_t)}\\|_2^2\\right]\\\\[.2cm]\n",
    "\\mbox{subject to} & x_{t+1} = Ax_t + B\\phi(x_t) + \\omega_t,\n",
    "\\end{array}\n",
    "\\label{eq:adp}\n",
    "\\end{equation}\n",
    "where $x_t\\in\\mathbf{R}^n$ is the state, $\\phi:\\mathbf{R}^n \\to \\mathcal U \\subseteq \\mathbf{R}^m$ is\n",
    "the policy, $\\mathcal U$ is a convex set representing the allowed set of controls,\n",
    "and $\\omega_t\\in\\Omega$ is a (random, i.i.d.) disturbance.\n",
    "Here the variable is the policy $\\phi$, and the expectation is taken over\n",
    "disturbances and the initial state $x_0$. If $\\mathcal U$ is not an affine\n",
    "set, then this problem is in general very difficult to solve.\n",
    "\n",
    "A common heuristic for solving stochastic control problems is\n",
    "approximate dynamic programming (ADP), which parametrizes $\\phi$\n",
    "and replaces the minimization over functions $\\phi$ with a minimization over parameters.\n",
    "In this example, we take $\\mathcal U$ to be the unit ball and we represent $\\phi$\n",
    "as a particular quadratic *control-Lyapunov* policy.\n",
    "Evaluating $\\phi$ corresponds to solving the SOCP\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & u^T P u + x_t^T Q u + q^T u \\\\\n",
    "\\mbox{subject to} & \\|{u}\\|_2 \\leq 1,\n",
    "\\end{array}\n",
    "\\label{eq:policy}\n",
    "\\end{equation}\n",
    "with variable $u$ and parameters $P$, $Q$, $q$, and $x_t$. We can run\n",
    "stochastic gradient descent (SGD) on $P$, $Q$, and $q$ to\n",
    "approximately solve the original problem, which requires requires differentiating\n",
    "through the quadratic policy. Note that if $u$ were unconstrained, the original problem\n",
    "could be solved exactly, via linear quadratic regulator (LQR) theory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "from scipy.linalg import solve_discrete_are\n",
    "from scipy.linalg import sqrtm\n",
    "\n",
    "\n",
    "from cvxpylayers.tensorflow.cvxpylayer import CvxpyLayer\n",
    "\n",
    "\n",
    "# Generate data\n",
    "tf.random.set_seed(1)\n",
    "np.random.seed(1)\n",
    "\n",
    "n = 2\n",
    "m = 3\n",
    "\n",
    "A = np.eye(n) + 1e-2 * np.random.randn(n, n)\n",
    "B = 1e-2 / 3 * np.random.randn(n, m)\n",
    "Q = np.eye(n)\n",
    "R = np.eye(m)\n",
    "\n",
    "# Compute LQR control policy\n",
    "P_lqr = solve_discrete_are(A, B, Q, R)\n",
    "P = R + B.T@P_lqr@B\n",
    "P_sqrt_lqr = sqrtm(P)\n",
    "\n",
    "# Construct CVXPY problem and layer\n",
    "x_cvxpy = cp.Parameter((n, 1))\n",
    "P_sqrt_cvxpy = cp.Parameter((m, m))\n",
    "P_21_cvxpy = cp.Parameter((n, m))\n",
    "q_cvxpy = cp.Parameter((m, 1))\n",
    "\n",
    "u_cvxpy = cp.Variable((m, 1))\n",
    "y_cvxpy = cp.Variable((n, 1))\n",
    "\n",
    "objective = .5 * cp.sum_squares(P_sqrt_cvxpy @ u_cvxpy) + x_cvxpy.T @ y_cvxpy + q_cvxpy.T @ u_cvxpy\n",
    "problem = cp.Problem(cp.Minimize(objective), [cp.norm(u_cvxpy) <= 1, y_cvxpy == P_21_cvxpy @ u_cvxpy])\n",
    "assert problem.is_dpp()\n",
    "policy = CvxpyLayer(\n",
    "    problem, [x_cvxpy, P_sqrt_cvxpy, P_21_cvxpy, q_cvxpy], [u_cvxpy])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we train the policy and plot the estimated average cost for each iteration of\n",
    "SGD for a numerical example, with $x \\in\n",
    "\\mathbf{R}^2$ and $u \\in \\mathbf{R}^3$, a time horizon of $T=25$, and a batch\n",
    "size of $8$. We initialize our policy's parameters with the LQR solution,\n",
    "ignoring the constraint on $u$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(iter 0) loss: 1.18634 \n",
      "(iter 1) loss: 1.17733 \n",
      "(iter 2) loss: 1.16115 \n",
      "(iter 3) loss: 1.13998 \n",
      "(iter 4) loss: 1.11706 \n",
      "(iter 5) loss: 1.09234 \n",
      "(iter 6) loss: 1.06502 \n",
      "(iter 7) loss: 1.03633 \n",
      "(iter 8) loss: 1.00895 \n",
      "(iter 9) loss: 0.982198 \n",
      "(iter 10) loss: 0.954998 \n",
      "(iter 11) loss: 0.926645 \n",
      "(iter 12) loss: 0.900704 \n",
      "(iter 13) loss: 0.874988 \n",
      "(iter 14) loss: 0.851813 \n",
      "(iter 15) loss: 0.829279 \n",
      "(iter 16) loss: 0.809518 \n",
      "(iter 17) loss: 0.79314 \n",
      "(iter 18) loss: 0.77915 \n",
      "(iter 19) loss: 0.767274 \n",
      "(iter 20) loss: 0.756987 \n",
      "(iter 21) loss: 0.74633 \n",
      "(iter 22) loss: 0.735943 \n",
      "(iter 23) loss: 0.725892 \n",
      "(iter 24) loss: 0.716452 \n",
      "(iter 25) loss: 0.70797 \n",
      "(iter 26) loss: 0.700616 \n",
      "(iter 27) loss: 0.694235 \n",
      "(iter 28) loss: 0.688617 \n",
      "(iter 29) loss: 0.683476 \n",
      "(iter 30) loss: 0.678953 \n",
      "(iter 31) loss: 0.67472 \n",
      "(iter 32) loss: 0.670976 \n",
      "(iter 33) loss: 0.667654 \n",
      "(iter 34) loss: 0.664699 \n",
      "(iter 35) loss: 0.662062 \n",
      "(iter 36) loss: 0.659702 \n",
      "(iter 37) loss: 0.657582 \n",
      "(iter 38) loss: 0.655674 \n",
      "(iter 39) loss: 0.653949 \n",
      "(iter 40) loss: 0.652387 \n",
      "(iter 41) loss: 0.650967 \n",
      "(iter 42) loss: 0.649673 \n",
      "(iter 43) loss: 0.648491 \n",
      "(iter 44) loss: 0.647407 \n",
      "(iter 45) loss: 0.646412 \n",
      "(iter 46) loss: 0.645496 \n",
      "(iter 47) loss: 0.64465 \n",
      "(iter 48) loss: 0.643867 \n",
      "(iter 49) loss: 0.64314 \n",
      "(iter 50) loss: 0.642464 \n",
      "(iter 51) loss: 0.641835 \n",
      "(iter 52) loss: 0.641246 \n",
      "(iter 53) loss: 0.640696 \n",
      "(iter 54) loss: 0.640179 \n",
      "(iter 55) loss: 0.639693 \n",
      "(iter 56) loss: 0.639236 \n",
      "(iter 57) loss: 0.638805 \n",
      "(iter 58) loss: 0.638397 \n",
      "(iter 59) loss: 0.63801 \n",
      "(iter 60) loss: 0.637644 \n",
      "(iter 61) loss: 0.637296 \n",
      "(iter 62) loss: 0.636964 \n",
      "(iter 63) loss: 0.636648 \n",
      "(iter 64) loss: 0.636346 \n",
      "(iter 65) loss: 0.636058 \n",
      "(iter 66) loss: 0.635782 \n",
      "(iter 67) loss: 0.635517 \n",
      "(iter 68) loss: 0.635262 \n",
      "(iter 69) loss: 0.635018 \n",
      "(iter 70) loss: 0.634782 \n",
      "(iter 71) loss: 0.634555 \n",
      "(iter 72) loss: 0.634336 \n",
      "(iter 73) loss: 0.634125 \n",
      "(iter 74) loss: 0.63392 \n",
      "(iter 75) loss: 0.633722 \n",
      "(iter 76) loss: 0.633531 \n",
      "(iter 77) loss: 0.633345 \n",
      "(iter 78) loss: 0.633164 \n",
      "(iter 79) loss: 0.632989 \n",
      "(iter 80) loss: 0.632819 \n",
      "(iter 81) loss: 0.632653 \n",
      "(iter 82) loss: 0.632492 \n",
      "(iter 83) loss: 0.632335 \n",
      "(iter 84) loss: 0.632182 \n",
      "(iter 85) loss: 0.632033 \n",
      "(iter 86) loss: 0.631887 \n",
      "(iter 87) loss: 0.631745 \n",
      "(iter 88) loss: 0.631606 \n",
      "(iter 89) loss: 0.631471 \n",
      "(iter 90) loss: 0.631339 \n",
      "(iter 91) loss: 0.631209 \n",
      "(iter 92) loss: 0.631082 \n",
      "(iter 93) loss: 0.630958 \n",
      "(iter 94) loss: 0.630837 \n",
      "(iter 95) loss: 0.630718 \n",
      "(iter 96) loss: 0.630601 \n",
      "(iter 97) loss: 0.630487 \n",
      "(iter 98) loss: 0.630375 \n",
      "(iter 99) loss: 0.630265 \n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def train(iters):\n",
    "    # Initialize with LQR control lyapunov function\n",
    "    P_sqrt = tf.Variable(P_sqrt_lqr)\n",
    "    P_21 = tf.Variable(A.T @ P_lqr @ B)\n",
    "    q = tf.Variable(tf.zeros((m, 1), dtype=tf.float64))\n",
    "    variables = [P_sqrt, P_21, q]\n",
    "    A_tf, B_tf, Q_tf, R_tf = map(tf.constant, [A, B, Q, R])\n",
    "\n",
    "    def g(x, u):\n",
    "        return tf.squeeze(tf.transpose(x) @ Q_tf @ x + tf.transpose(u) @ R_tf @ u)\n",
    "\n",
    "    def evaluate(x0, P_sqrt, P_21, q, T):\n",
    "        x = x0\n",
    "        cost = 0.\n",
    "        for _ in range(T):\n",
    "            u, = policy(x, P_sqrt, P_21, q)\n",
    "            cost += g(x, u) / T\n",
    "            x = A_tf @ x + B_tf @ u + .2 * tf.random.normal((n, 1), dtype=tf.float64)\n",
    "        return cost\n",
    "\n",
    "    def eval_loss(N=8, T=25):\n",
    "        return sum([evaluate(tf.zeros((n, 1), dtype=tf.float64), P_sqrt, P_21, q, T=T)\n",
    "                    for _ in range(N)]) / N\n",
    "\n",
    "    results = []\n",
    "    optimizer = tf.keras.optimizers.SGD(learning_rate=0.02, momentum=0.9)\n",
    "    for i in range(iters):\n",
    "        tf.random.set_seed(1)\n",
    "        np.random.seed(1)\n",
    "        with tf.GradientTape() as tape:\n",
    "            loss = eval_loss()\n",
    "        gradients = tape.gradient(loss, variables)\n",
    "        optimizer.apply_gradients(zip(gradients, variables))\n",
    "        results.append(loss.numpy())\n",
    "        print(\"(iter %d) loss: %g \" % (i, results[-1]))\n",
    "    return results\n",
    "\n",
    "\n",
    "results = train(iters=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhU5Zn38e9d1RtN0zTQG9BsCgKCKIjEPRo31ERNNIlGk8zExDiTbTLJTJJ5M76JM85Mksk6Y5JXjUs2E2NM3LegCa7RBhFBQBAUWpZuduiml6q63z/qtLbYDQX06dNd5/e5rrqqzlJV9/Fg/fo855znMXdHRETiKxF1ASIiEi0FgYhIzCkIRERiTkEgIhJzCgIRkZgriLqAA1VZWenjx4+PugwRkQFlwYIFm929qrtlAy4Ixo8fT319fdRliIgMKGb2ek/LQmsaMrObzazRzJb0sPxyM1scPJ42s6PDqkVERHoW5jmCW4G5+1i+Bni3u88A/g24IcRaRESkB6E1Dbn7fDMbv4/lT3eZfBaoC6sWERHpWX+5auhK4MGeFprZVWZWb2b1TU1NfViWiEj+izwIzOx0skHwlZ7Wcfcb3H22u8+uqur2pLeIiBykSK8aMrMZwE3Aue6+JcpaRETiKrIjAjMbC9wFfNTdX4mqDhGRuAvz8tHbgWeAyWbWYGZXmtnVZnZ1sMo1wAjgx2a2yMxCvTlgy+42vnnvUtpS6TC/RkRkwAnzqqHL9rP8k8Anw/r+vT27eiu3PPUaDdv28OPLZ1GYjPz0iIhIvxCbX8PzZ4zkmxdM49GXN/HF3y4indGAPCIiMAC7mDgUHz9xPK0daf7zweWUFCb59sUzSCQs6rJERCIVqyAA+PS7D6elPc0P563k+MNGcMmxuo9NROItNk1DXX3hjEkcNXoo33/0FZ08FpHYi2UQJBLGV+ZO4Y3te/jls2ujLkdEJFKxDAKAkydVcvLESq5/fBW7WjuiLkdEJDKxDQKAr8ydwtbmdm6cvzrqUkREIhPrIDiqbijnzxjJTU+uoWlXW9TliIhEItZBAPCls45gT0eaXz7b4+A9IiJ5LfZBcFhVGadOquI3z6+lI52JuhwRkT4X+yAAuOL4cWza2ca8ZY1RlyIi0ucUBMDpk6sYObSEX/1VzUMiEj8KAqAgmeCyOWN5YuVmXtvcHHU5IiJ9SkEQ+PBxY0gmjF8/pxvMRCReFASBmvISzj6yht/Vr6O1Q91OiEh8KAi6uPxd49jW0sFDSzZGXYqISJ9REHRx4uEjGDN8EL9bsC7qUkRE+oyCoItEwrh4Vh1Pv7qFhm0tUZcjItInFAR7uXhWHe5w18I3oi5FRKRPKAj2MmZ4KScePoI7FzSQ0XCWIhIDCoJuXHJsHWu3tvDca1ujLkVEJHQKgm6cO30kZcUF/K6+IepSRERCpyDoxqCiJO+dMZIHl2yguS0VdTkiIqFSEPTgg7PraGlPc/9LG6IuRUQkVAqCHswaO4xxI0q5Z9H6qEsREQmVgqAHZsYFR4/i6Vc307irNepyRERCoyDYhwuPGUXG4f7Fah4SkfwVWhCY2c1m1mhmS3pYPsXMnjGzNjP7clh1HIqJ1UOYOrKce15U85CI5K8wjwhuBebuY/lW4PPAf4dYwyG78JhRvLB2O2u3qMsJEclPoQWBu88n+2Pf0/JGd38e6Airht7wvqNHAXDvYh0ViEh+GhDnCMzsKjOrN7P6pqamPv3u0RWDOG78MO5epL6HRCQ/DYggcPcb3H22u8+uqqrq8++/4OhRvLJpN8s37uzz7xYRCduACIKonXfUSJIJ0z0FIpKXFAQ5GFFWzPGHDefhpRq5TETyT5iXj94OPANMNrMGM7vSzK42s6uD5bVm1gD8I/D1YJ3ysOo5VOdMq+XVpmZWNe6KuhQRkV5VENYHu/tl+1m+EagL6/t729lH1nLN3Ut5eOkmJlYPibocEZFeo6ahHNUOLeHoMRVqHhKRvKMgOADnTKthccMO1m/fE3UpIiK9RkFwAM6ZVgvAIzoqEJE8oiA4AIdXlTGxuoyHl26KuhQRkV6jIDhAc6fV8txrW9nW3B51KSIivUJBcIDOmVZLOuP8aZmOCkQkPygIDtD00eWMGlrCIy8rCEQkPygIDpCZccbUGp5cuZnWjnTU5YiIHDIFwUE4Y2o1ezrSPPPqlqhLERE5ZAqCg3DC4SMYXJTkUZ0nEJE8oCA4CMUFSU6ZVMVjyxpx96jLERE5JAqCg3TG1Go27mxl6XqNUSAiA5uC4CCdPqUaM3QZqYgMeAqCg1RZVszMMRXMW9YYdSkiIodEQXAIzphaw0tv7GDjjtaoSxEROWgKgkNw5tQaAOYtV/OQiAxcCoJDcERNGWOGD+JPustYRAYwBcEhMDPOmlrLU6u2sLstFXU5IiIHRUFwiM6ZVkN7OsOfV+iksYgMTAqCQzR7/HBGDC7SGAUiMmApCA5RMmGcObWGx5c30pZSJ3QiMvAoCHrB2dNq2N2WUid0IjIgKQh6wUkTKxlclFTzkIgMSAqCXlBSmOS0ydU8+vIm0hl1QiciA4uCoJecPa2GzbvbeGHttqhLERE5IAqCXnL6lGoKk8bDSzdGXYqIyAEJLQjM7GYzazSzJT0sNzP7kZmtMrPFZjYrrFr6QnlJISceXsnDSzdpjAIRGVDCPCK4FZi7j+XnApOCx1XAT0KspU/MnV7L2q0tvLxBYxSIyMARWhC4+3xg6z5WuRD4uWc9C1SY2ciw6ukLZx9ZQ8LgoSVqHhKRgSPKcwSjgXVdphuCee9gZleZWb2Z1Tc1NfVJcQdjRFkxcyYMVxCIyIASZRBYN/O6bVx39xvcfba7z66qqgq5rENz7vSRrGzczarG3VGXIiKSk/0GgZl9MJd5B6EBGNNlug5Y3wufG6lzptUC8NCSDRFXIiKSm1yOCL6W47wDdQ/wseDqoeOBHe4+4H89a4eWMGtsBQ+qeUhEBoiCnhaY2bnAecBoM/tRl0XlwH473zez24HTgEozawD+L1AI4O4/BR4IPn8V0AL87cFtQv9z7vSRXPfAMtZuaWHsiNKoyxER2aceg4BsM009cAGwoMv8XcAX9/fB7n7ZfpY78Jkcahxw5k6v5boHlvHw0o186tTDoi5HRGSfegwCd38ReNHMfu3uHQBmNgwY4+7qR2EfxgwvZdqoch5cskFBICL9Xi7nCB41s3IzGw68CNxiZt8Lua4B79zptSxcu51NO1ujLkVEZJ9yCYKh7r4T+ABwi7sfC5wZblkDX+fVQ49oYHsR6edyCYKC4I7fDwH3hVxP3phYXcZhlYN5RJ3QiUg/l0sQXAs8DLzq7s+b2WHAynDLGvjMjHOm1/LMq1vY0dIRdTkiIj3abxC4++/cfYa7/10wvdrdLw6/tIHvnGm1pDLOvOVqHhKR/iuXO4vrzOwPQZfSm8zs92ZW1xfFDXQzRg+ltrxEfQ+JSL+WS9PQLWTvAh5FtlO4e4N5sh+JhHHOtBrmr2yipX2/9+CJiEQilyCocvdb3D0VPG4F+nfPb/3IOdNqae3IMP+V/ttrqojEWy5BsNnMrjCzZPC4AtgSdmH5Ys6E4VSUFvLwUp0nEJH+KZcg+ATZS0c3AhuAS4J5koOCZIIzp9bwp2Wb6Ehnoi5HROQdcrlqaK27X+DuVe5e7e4XufvrfVFcvjjryBp2taZ4/rV9DdgmIhKNXK4aus3MKrpMDzOzm8MtK7+cMqmSooIEf3q5MepSRETeIZemoRnuvr1zIuhwbmZ4JeWf0qICTp5YyaPLNpLtdFVEpP/IJQgSQa+jAASdz+2r+2rpxplTa1i3dQ8rNYSliPQzuQTBd4GnzezfzOxa4Gng2+GWlX/OmFoNwKPqhE5E+plcThb/HLgY2AQ0AR9w91+EXVi+qSkv4ei6oQoCEel3cmricfeXgZdDriXvnTm1hu8++gqNu1qpHlISdTkiIkBuTUPSS86aVgPAY8t09ZCI9B8Kgj40uWYIdcMG8adlah4Skf4jpyAws3FmdmbwepCZDQm3rPxkZpw5tYYnVm6muU2d0IlI/5DLDWWfAu4E/l8wqw74Y5hF5bNzp9fSlsrw2HI1D4lI/5DLEcFngJOAnQDuvhKoDrOofDZ7/HCqhhTzwEsboi5FRATILQja3L29c8LMCgDdHnuQkgnj3Om1PL6iUc1DItIv5BIEfzGzfwEGmdlZwO/IDk4jB+n8o0bS2qHmIRHpH3IJgq+SvZHsJeDTwAPA18MsKt+peUhE+pP93lDm7hngxuAhvaCzeei3z6+juS3F4GJ13SQi0cnlqqGXzGzxXo8nzOz7ZjZiP++da2YrzGyVmX21m+XjzGxe8Jl/NrO6Q9mYgeS8o0bq6iER6RdyaRp6ELgfuDx43AvMJzti2a09vcnMksD1wLnAkcBlZnbkXqv9N/Bzd58BXAv85wHWP2AdFzQP3b9YzUMiEq1c2iROcveTuky/ZGZPuftJwfjFPZkDrHL31QBm9hvgQt7eZ9GRwBeD148To/sT1DwkIv1FLkcEZWb2rs4JM5sDlAWT+7r+cTSwrst0QzCvqxfJ9mwK8H5gSHfNTWZ2lZnVm1l9U1NTDiUPDOcHzUPz1DwkIhHKJQg+CdxkZmvM7DXgJuBTZjaYfTflWDfz9r7/4MvAu83sBeDdwBt0Ey7ufoO7z3b32VVVVTmUPDDMHj+c6iHF3L94fdSliEiM5XLV0PPAUWY2FLCuw1YCd+zjrQ3AmC7TdcDbfvHcfT3wAQAzKwMudvcdOdY+4CUTxnlHjeTXz61ld1uKMjUPiUgEcu107nyy9xB83syuMbNrcnjb88AkM5tgZkXApcA9e31upZl11vA14ObcS88P588YSXsqwzz1SCoiEcnl8tGfAh8GPke2ueeDwLj9vc/dU8BngYeBZcAd7r7UzK41swuC1U4DVpjZK0ANcN3BbMRAduzYYdSWl3Cfrh4SkYjk0hZxorvPMLPF7v5NM/sucFcuH+7uD5C9E7nrvGu6vL6TbM+msZUImod++ezr7GrtYEhJYdQliUjM5NI01Bo8t5jZKKADmBBeSfFz/oyRtKczGs9YRCKRSxDca2YVwHeAhcBrwO1hFhU3M8dUMGpoiW4uE5FI7LNpKDiROy+4Uuj3ZnYfUBKnK3v6Qmfz0G3PvMaOlg6Glqp5SET6zj6PCIIO577bZbpNIRCOC48ZTUfaue8l3VMgIn0rl6ahR8zsYjPr7gYx6SXTR5dzRE0Zv1/QEHUpIhIzuQTBP5IdjKbdzHaa2S4z2xlyXbFjZlxybB0L127n1abdUZcjIjGy3yBw9yHunnD3QncvD6bL+6K4uLnomNEkDO5aqKMCEek7udxQZmZ2hZn9azA9Juh4TnpZdXkJpx5RxV0L3yCd0bDQItI3cmka+jFwAvCRYHo32XEGJASXHFvHhh2tPPPqlqhLEZGYyCUI3uXunyG4sczdtwFFoVYVY2dOraG8pIA7F6zb/8oiIr0glyDoCEYbcwAzqwIyoVYVYyWFSd539CgeWrqRXa0dUZcjIjGQSxD8CPgDUG1m1wFPAv8RalUxd8mxdbR2ZHSnsYj0iVzGI/iVmS0AziDb++hF7r4s9Mpi7JgxFUyqLuOO+nVcOmds1OWISJ7L5aqhHwLD3f16d/9fhUD4zIwPzR7DwrXbWdW4K+pyRCTP5dI0tBD4upmtMrPvmNnssIsSeP+s0RQkjDvqdU+BiIQrlxvKbnP384A5wCvAt8xsZeiVxVxlWTFnTK3mroUNdKR1bl5EwpPTUJWBicAUYDywPJRq5G0+NHsMm3e38/jyxqhLEZE8lss5gs4jgGuBpcCx7v6+0CsT3n1EFVVDitU8JCKhymWoyjXACe6+Oexi5O0KkgkunlXHjU+spnFnK9XlJVGXJCJ5KJdzBD8F0mY2x8xO7Xz0QW0CfPi4MaQzzu3P6U5jEQlHLk1DnwTmAw8D3wyevxFuWdJpQuVgTp9cxS+efZ22VDrqckQkD+VysvgLwHHA6+5+OjATaAq1KnmbT5w8gc2727j3Rd1pLCK9L5cgaHX3VgAzK3b35cDkcMuSrk6eWMkRNWXc/OQa3NU9tYj0rlyCoMHMKoA/Ao+a2d2ABtbtQ2bGJ06awMsbdvLXNVujLkdE8kwuJ4vf7+7b3f0bwL8CPwMuCrswebuLZo5mWGkhNz+5JupSRCTPHMgNZbj7X9z9HndvD6sg6V5JYZLL3zWOR5dt4vUtzVGXIyJ55ICC4ECZ2VwzWxH0U/TVbpaPNbPHzewFM1tsZueFWc9A99ETxlGYTPC/j62KuhQRySOhBUEwmM31wLnAkcBlZnbkXqt9HbjD3WcCl5IdFlN6UFNewkePH8fvFzaoV1IR6TVhHhHMAVa5++qgKek3wIV7reNAefB6KDoJvV9/f9rhDCpM8t1HXom6FBHJE2EGwWig6+2wDcG8rr4BXGFmDcADwOdCrCcvjCgr5pOnHMaDSzayuGF71OWISB4IMwism3l7XwR/GXCru9cB5wG/MLN31GRmV5lZvZnVNzXpXrZPnjKB4YOL+M7DK6IuRUTyQJhB0ACM6TJdxzubfq4E7gBw92eAEqBy7w9y9xvcfba7z66qqgqp3IFjSEkhf3/a4TyxcjNPrVJfgCJyaMIMgueBSWY2wcyKyJ4MvmevddaSHQsZM5tKNgj0J38Orjh+HGOGD+Kau5fQntLANSJy8EILAndPAZ8l20ndMrJXBy01s2vN7IJgtS8BnzKzF4Hbgb9x9aGQk5LCJNdeMJ1Xm5q58YnVUZcjIgNYLuMRHDR3f4DsSeCu867p8vpl4KQwa8hnp0+p5ryjavnRvJW8b8Yoxo4ojbokERmAQr2hTMJ3zXunUZAw/vXuJeqQTkQOioJggKsdWsKXzp7MX15p4t7F6qZaRA6cgiAPfOyEcRwzpoKv/+ElGra1RF2OiAwwCoI8UJBM8KNLZ5Jx+IffLCKV1lVEIpI7BUGeGDuilOveP53617fxw3kroy5HRAYQBUEeufCY0Xzw2Dr+9/FVPK0bzUQkRwqCPPONC6ZxWOVgPvPrhRq3QERyoiDIM4OLC7jp48eRcbjytnp2tnZEXZKI9HMKgjw0oXIwP73iWF7b3MxnfrVQJ49FZJ8UBHnqhMNHcN37p/PEys38+/3Loi5HRPqxULuYkGh9+LixvLJpNz97cg1HjiznQ8eN2f+bRCR2dESQ57527hROmVTJ1/+4hAWvb4u6HBHphxQEea4gmeB/LpvJyIoSrv7lAjbuaI26JBHpZxQEMVBRWsSNH5tNS1uKT/+intaOdNQliUg/oiCIiSNqhvCDS2ey+I0d/NOdi9VTqYi8SUEQI2cdWcM/nzOFe19cz/88tirqckSkn9BVQzFz9bsPY2XjLr736CtMrC7jvKNGRl2SiERMRwQxY2b8x/uPYtbYCr7420W8sFZXEonEnYIghkoKk9zwsdlUlxfzqZ/Xs26rxjAQiTMFQUxVlhVzy98cR3sqw9/e+jw79qhPIpG4UhDE2MTqIfz0o8fy+pZmrv7FAtpSuqxUJI4UBDF34uGVfOviGTyzegtf/O0i0hldVioSN7pqSPjArDq27G7nugeWMax0Cf9+0XTMLOqyRKSPKAgEgE+dehhbmtv56V9eZfjgIr509uSoSxKRPqIgkDd9Ze5ktjW38z+PrSJhxj+cOUlHBiIxoCCQN5kZ//GBo8i488N5K0llMnz57MkKA5E8pyCQt0kmjG9dPIOCZILrH3+V9lSGfzlvqsJAJI+FGgRmNhf4IZAEbnL3/9pr+feB04PJUqDa3SvCrEn2L5EwrrtoOoVJ48Yn1rBpZxvfvmQGJYXJqEsTkRCEFgRmlgSuB84CGoDnzewed3+5cx13/2KX9T8HzAyrHjkwiYTxzQumUTu0hG8/tIKGbS3c+LHZjCgrjro0EellYd5HMAdY5e6r3b0d+A1w4T7Wvwy4PcR65ACZGX9/2kR+fPkslq7fyYXXP8Xihu1RlyUivSzMIBgNrOsy3RDMewczGwdMAB7rYflVZlZvZvVNTU29Xqjs23lHjeSOT59AJuNc/JOnuemJ1RrPQCSPhBkE3Z1d7OnX41LgTnfvto8Dd7/B3We7++yqqqpeK1Byd/SYCh74wimcNrmaf79/GVfeVk/jLg17KZIPwgyCBmBMl+k6YH0P616KmoX6vYrSIm746LF8431H8uSqzZz1vfn8rn6djg5EBrgwg+B5YJKZTTCzIrI/9vfsvZKZTQaGAc+EWIv0EjPjb06awINfOIUjasr4pzsX87Gbn+O1zc1RlyYiBym0IHD3FPBZ4GFgGXCHuy81s2vN7IIuq14G/Mb1Z+WAcnhVGb+96gSuvXAaC1/fxtnfn893H1nBnnb1YCoy0NhA+/2dPXu219fXR12GdLFpZyv/+cAy/rhoPaMrBvGls4/gwmNGk0zoJjSR/sLMFrj77O6WqRtqOWQ15SX84NKZ3PHpE6goLeQf73iRuT+Yz0NLNur8gcgAoCCQXjNnwnDu/ezJ/PjyWWTcufqXCzj/R0/y4EsbyGicA5F+S01DEopUOsPdi9Zz/eOrWL25mSNqyrjq1MN539EjKS5QVxUifW1fTUMKAglVOuPct3g9P378VVZs2kVlWREfedc4PjJnLLVDS6IuTyQ2FAQSOXfnqVVbuOWpNTy2ohED3n1EFR+aPYYzptZQVKBWSpEwKQikX3l9SzN31K/jzgUNbNrZxtBBhZwzrYb3zhjFiYePoCCpUBDpbQoC6ZfSGWf+yibuXbSeR17exO62FEMHFfKeKdWcObWGU46opLykMOoyRfLCvoJAA9NIZJIJ4/TJ1Zw+uZrWjjR/XtHEIy9v5PHljfzhhTdIJoyZYyo49YgqTpo4gqNGV6gJSSQEOiKQfiedcRa8vo2/vNLIEys389IbO3CHQYVJZo2rYM74EcwaV8HRYyp0xCCSIzUNyYC2tbmdv67ewl/XbOXZ1VtYsWkX7mCW7epixuihTA8eU0YOUTiIdENBIHllZ2sHi9ft4IW121i0bjsvvbGDxl1tby6vGzaIKbVDmFQzhEnVZUysLuOwqjLKitUSKvGlcwSSV8pLCjl5UiUnT6p8c17jzlaWrN/Bsg27WL5xF8s37OQvrzTRkX7rD53qIcVMqBzMuBGljBsxmLHDSxk7vJS6YYMYPrgIM/WNJPGkIJC8UF1ewnvKS3jPlJo353WkM7y+pYVVjbtYvbmZNU3NrNnczOMrmmja1fC29w8qTDKyooTRFYMYObSE2qGDqC0vYeTQEqrLi6kaUsyIwcXqSE/ykoJA8lZhMsHEoGloby3tKdZubaFh6x4atrWwbtseNuzYwxvbW1m+sYnNu9vYu9U0mTCGDy6iqqyYyiHFVJYVMWJwESPKihk+uIjhpUUMG1zEsNJCKkqLGDqoUMEhA4KCQGKptKiAKbXlTKkt73Z5RzpD0642NuxopWlXK4272mjc2UbTrjY2784+Xm3czZbmNlo7Mt1+hhkMKS6gorSIitJChg4qpLykkPJBhZQPKqC8pJAhJQXZR3EhZSUFlBVnpwcXZ18XFyTUZCWhUxCIdKMwmWBUxSBGVQza77ot7Sm27G5nW0s7W5uzz9tbOtjW0sGOlnZ27Olg+54OduzpYP32PexsTbFjTwftqe4DpKuChFFalKSsOBsOpcUFDC5KUlpUQGlRktKiJIOC59KiAkoKg3mFSUoKs8uyrxNvzisuTFBSmKSkIElh0hQ0oiAQOVSlRQWUDi9gzPDSA3pfWyrNrtYUO/d0sLstxe7WFDtbUzS3pWhuT7Gr83Vbit1tafZ0pGhuS9PclmJ7yx5a2lM0t6fZ056mpT3FwfT0nTAoLsgGRUlhkuKCt56LChIUF3R9nX0uKkhQlMwGSlHy7csKk9l5hQWdy4yiZDZw3pqXXa8waRQlExQEr7PzEmpOi4CCQCQixQVJisuSVJYVH/JnuTvt6Qyt7RlaOlK0tKdp7cg+9rRn2NPx1nRrR5q2VCZ4naEtlX1u7UjTns7Q1pGhNZWmPZV937aWdtpTGdpSGdpTmWCd7Lpdr8rqLQmDgiBQCjoDImFvC4yCpFGQyE4XJDqn31qnc15hIkEyaRQmjGSwfrJzvYQFyxLBvGBZ4q33d04nE4ng2d56Tr59frLL8oR1/bwESct+V9LeWi9h9JujMQWBSB4ws+Cv9yRD6bsb6jIZpyMTBEQQEu2pDB3pt4IjlXE6Uhna0hk6UtnwaE+n6Ug7HcH6qXT2czpS2XkdmWBeEDbZ57fmpTL+5nQqk2FPh5POvLUsFbwvncku73zduTzdTwZKSibeGQ4FyUQ2SPYKmITBZXPG8slTDuv1OhQEInLQEgmjOJEccIMNuXeGRPaRDoLozXnp7OtseDgZD9YLQiXz5nszpDOQzrwVMJ2fken6Wc6b0xl3Umkn7W+tkw4+J/s9weuu6wef1xtHj91REIhI7FjQdDPA8is06spRRCTmFAQiIjGnIBARiTkFgYhIzCkIRERiTkEgIhJzCgIRkZhTEIiIxNyAG6rSzJqA1w/y7ZXA5l4sZ6CI43bHcZshntsdx22GA9/uce5e1d2CARcEh8LM6nsaszOfxXG747jNEM/tjuM2Q+9ut5qGRERiTkEgIhJzcQuCG6IuICJx3O44bjPEc7vjuM3Qi9sdq3MEIiLyTnE7IhARkb0oCEREYi42QWBmc81shZmtMrOvRl1PGMxsjJk9bmbLzGypmX0hmD/czB41s5XB87Coaw2DmSXN7AUzuy+YnmBmfw22+7dmVhR1jb3JzCrM7E4zWx7s8xPisK/N7IvBv+8lZna7mZXk4742s5vNrNHMlnSZ1+3+tawfBb9vi81s1oF8VyyCwMySwPXAucCRwGVmdmS0VYUiBXzJ3acCxwOfCbbzq8A8d58EzAum89EXgGVdpr8FfD/Y7m3AlZFUFZ4fAg+5+xTgaLLbntf72sxGA58HZrv7dCAJXEp+7utbgbl7zetp/54LTAoeVwE/OZAvikUQAHOAVe6+2t3bgd8AF0ZcU69z9w3uvjB4vYvsD8Nostt6W7DabcBF0VQYHjOrA84HbgqmDXgPcGewSuPJapUAAAQ/SURBVF5tt5mVA6cCPwNw93Z3304M9jXZIXYHmVkBUApsIA/3tbvPB7buNbun/Xsh8HPPehaoMLORuX5XXIJgNLCuy3RDMC9vmdl4YCbwV6DG3TdANiyA6ugqC80PgH8GMsH0CGC7u6eC6Xzb54cBTcAtQXPYTWY2mDzf1+7+BvDfwFqyAbADWEB+7+uuetq/h/QbF5cgsG7m5e11s2ZWBvwe+Ad33xl1PWEzs/cCje6+oOvsblbNp31eAMwCfuLuM4Fm8qwZqDtBm/iFwARgFDCYbLPI3vJpX+fikP69xyUIGoAxXabrgPUR1RIqMyskGwK/cve7gtmbOg8Tg+fGqOoLyUnABWb2Gtlmv/eQPUKoCJoPIP/2eQPQ4O5/DabvJBsM+b6vzwTWuHuTu3cAdwEnkt/7uque9u8h/cbFJQieByYFVxYUkT25dE/ENfW6oF38Z8Ayd/9el0X3AB8PXn8cuLuvawuTu3/N3evcfTzZffuYu18OPA5cEqyWV9vt7huBdWY2OZh1BvAyeb6vyTYJHW9mpcG/987tztt9vZee9u89wMeCq4eOB3Z0NiHlxN1j8QDOA14BXgX+T9T1hLSNJ5M9HFwMLAoe55FtL58HrAyeh0dda4j/DU4D7gteHwY8B6wCfgcUR11fL2/rMUB9sL//CAyLw74GvgksB5YAvwCK83FfA7eTPQ/SQfYv/it72r9km4auD37fXiJ7VVXO36UuJkREYi4uTUMiItIDBYGISMwpCEREYk5BICIScwoCEZGYUxBIbJnZ08HzeDP7SC9/9r90910i/ZEuH5XYM7PTgC+7+3sP4D1Jd0/vY/ludy/rjfpEwqYjAoktM9sdvPwv4BQzWxT0dZ80s++Y2fNB3+6fDtY/LRjv4ddkb9rBzP5oZguC/vGvCub9F9neMReZ2a+6fldw5+d3gr70XzKzD3f57D93GV/gV8GdsyKhK9j/KiJ576t0OSIIftB3uPtxZlYMPGVmjwTrzgGmu/uaYPoT7r7VzAYBz5vZ7939q2b2WXc/ppvv+gDZO4KPBiqD98wPls0EppHtI+Ypsn0oPdn7myvydjoiEHmns8n227KIbDfeI8gO+AHwXJcQAPi8mb0IPEu2069J7NvJwO3unnb3TcBfgOO6fHaDu2fIdg8yvle2RmQ/dEQg8k4GfM7dH37bzOy5hOa9ps8ETnD3FjP7M1CSw2f3pK3L6zT6/1P6iI4IRGAXMKTL9MPA3wVdemNmRwSDvuxtKLAtCIEpZIcH7dTR+f69zAc+HJyHqCI7ythzvbIVIgdJf3GIZHvvTAVNPLeSHQt4PLAwOGHbRPdDHz4EXG1mi4EVZJuHOt0ALDazhZ7tErvTH4ATgBfJ9hT7z+6+MQgSkUjo8lERkZhT05CISMwpCEREYk5BICIScwoCEZGYUxCIiMScgkBEJOYUBCIiMff/AW29QgDRj+ZAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(results)\n",
    "plt.xlabel('iteration')\n",
    "plt.ylabel('average cost')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
